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SUMMARY 


A code has been developed for simulating high Reynolds number transonic 
flow fields of arbitrary configuration. This code, in conjunction with labo- 
ratory experiments, is being used to devise and test turbulence transport 
models which may be suitable in the prediction of such flow fields, with par- 
ticular emphasis on regions of flow separation. The solutions describe the 
flow field, including both the shock- induced and tralling-edge separation 
regions, in sufficient detail to provide the profile and friction drag. 
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INTRODUCTION 


Transonic flow fields over airfoil configurations are characterized by 
regions of subsonic flow and regions of supersonic flow. These regions are 
often separated by standing shock waves. If these shocks are strong enough, 
flow separation will occur where the shock impinges on the airfoil surface. 

In addition, separation is likely to occur near the airfoil trailing edge. 

If the Reynolds number is large, the flow will be turbulent near the airfoil 
surface and, in all cases, will be turbulent in the far wake. 

To simulate such flow fields numerically, it is necessary to use a set of 
equations capable of supporting a description of all the above phenomena. In 
the present study the two-dimensional time-averaged Navier-Stokes equations 
for compressible flow are used. The Reynolds stresses are described by an 
algebraic eddy-viscosity model and the resulting system is solved by the 
second-order-accurate explicit difference method developed by MacCormack 
(refs. 1 and 2). 

This paper describes the progress made to date in the development of such 
a code, the procedure used to validate the code, and the adaptation of the 
code to advanced computers. Solutions are shown and compared with experiments 
for the flow field over an 18-percent-thick, biconvex, circular-arc airfoil 
at zero angle of attack for several values of free-stream Mach number and 
chord Reynolds number. 
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SYMBOLS 


c chord length 

Cf skin-friction coefficient 

Cp pressure coefficient 

e internal energy per unit volume 

k thermal conductivity coefficient 

Z mixing length 

lijo free-stream Mach number 

p pressure 

Pr Prandtl number 

Pr^- turbulent Prandtl number 

Re^ chord Reynolds number 

s scalar area 

S Integrated surface area 

t time 

T temperature 

u,v velocity vector components in x- and y-direction, respectively 

vol Integrated volume element 

p coefficient of molecular viscosity 

p mass density, 

a normal stress 

T shear stress 

Subscript: 

w wall surface 

Superscript ; 

+ denotes nondimenslonalized boundary-layer quantity 
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SIMULATION METHOD 


Governing Equations 


The flow field Is described by the two-dimensional compressible equations 
of motion for turbulent flow. Written In time-dependent Integral form, they 
are: 


where 
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and e , e are unit vectors, and n Is a unit normal vector. These equa- 
X y 

tlons can be solved In the orthogonal x,y coordinate system for an arbitrary 
quadrilateral volume element (sketch (a)) by application of the split Ly and 
Lx operators In the manner described by MacCormack (ref. 2). 
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Sketcli (a) — Quadrilateral volume element.' 


The Ly and Lx operators are. applied so that equation (1) is satisfied at 
each time step for each .cell of the nonofthogohal nonunifprm computational 
mesh. To evaluate the viscous derivatives for the nonorthogonal mesh, the' 
following transformation is appropriate (ref. 3): 

3(j) _ 3({) 34> 3r) 3({i _ 3(|) 3^ , 3(|> 3n 

3x ” 3g 3x 3n 3x 3y “ 3^ 3y 3n 3y 

where <{•. is a dummy dependent variable, and (C>n) are the local coordinates 
of the nonorthogonal mesh. In terms of the notation in sketch (b), 

... Ad)_Ay - A(}) Ay. A(J)_Ax - A* Ax^ 

3x ~ Ax_Ay - Ax Ay_ ^ 3y Ay_Ax - Ay Ax^ 

' ^ r) n . C w . • 5 n n C 
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Sketch (b) — Nonorthogonal mesh notation. 
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For evaluation on surfaces S3 and S4, the differences are defined as 


and 
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This treatment of the viscous 'deriyatiyes always results iii centered differ- - 
ences, '^fhtains' second-order accuracy,’ and 'provides consistent treatment of 
discontinuous boundary conditions (such as at the leading arid tfailirig edges-;' 
of the airfoil) . ■ ' ‘ • •.;'i /) 


Control Volume, Mesh, and Boundary Conditions 

the airfoil, initially at rest, is Impulsively started at- time zero at 
the desired free-stream Mach number and pressure. Figure 1 shows a'typical- ' 
control volume for which the flow-field development is followed in time. At 
a sufficient distance upstream of the leading edge (in this case six chord 
lengths), the flow is assumed uniform at the -free-stream conditions (U = Uoo) 
as it is along the far transverse boundary (again, six chord lengths away) . 

The downstream boundary is positioned far enough downstream of the trailing 
edge (nine chord lengths) so that all gradients in the flow direction may be 
assumed negligible (9u/8x = 0). The surface of the airfoil is impermeable, 
and "no slip" boundary conditions are assumed (u = v = 0) . The airfoil is 
assumed adiabatic (VT ' n = 0) , and the normal surface pressure gradient is 
zero (3p/8n = 0). Ahead of and behind the airfoil, the flow is symmetric. If 
the airfoil is thick, and the flow field is transonic, significant boundary- 
layer separation is likely. To simulate this phenomenon reliably for turbu- 
lent flow it is necessary to resolve the bouridary lay er to the sublayer scale. 
This sublayer scale is nearly proportional to 1//Re[. so that, for the high 
Reynolds number flows of Interest, the mesh resolution near the surface must 
be extremely f ine. As a rule of thumb, a first mesh spacing of 
Ayj^in = 2/3(c//Re^,) is adequate. 

The mesh used in the studies to-date is a 50x38 mesh. In the x-dlrectlon, 
the mesh is uniformly distributed over the surface of the airfoil (20 points) 
and is exponentially stretched ahead of (10 points) and behind the airfoil 
(20 points) . In the y-direction, a coarse mesh of 26 points is exponentially 
stretched away from the airfoil. The innermost region is further subdivided 
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Into a medium mesh of 10 exponentially stretched points and a -fine mesh of 
4 uniformly spaced points. 


. . Turbulence Model 

. The turbulence modeling Is Incorporated In the shear stress term 
and Tvv ..Intthe form of an eddyTvlscoslty. coefficient e as , , 
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Two different algebraic" models. 'have 'been considered » each expected to perform 
adequately ahead of the shock-wave — boundary-layer interaction regions but to 
exhibit different behavior when the flow separates. 

,v . , J..- ... .■ 

The first model considered (model' 1) Is a van Driest formulation (ref. 4) 
for the wall region and a plane mixing .formulation for the outer region. No 
special consideration is given to. separated' regions and the airfoil wake Is 
described by a Clauser (ref. 5) wake formulation. 
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Z = 0.07(6 - y^) 

where 6 Is the boundary-layer thickness and y^ Is the furthest point 
across the boundary layer where the velocity Is zero. (For attached boundary 
layers y<, = 0.) 

Airfoil Wake; 

e = 0. 001176p(6 - yQ)|“j “ 

where u^ Is the velocity at the edge of the wake and u^ is the velocity at 
the center line. 
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Prandtl Number 


Pr = Pr^ = 0.90 


The second model (model 2) uses a van Driest formulation for the wall 
region and a Clauser formulation with an Intermittency factor in the outer 
region. A limiting mlniinum^yalue; is "Mpo'sed. on^veloclty gradient used in the 
separation-bubble description and’ the^ airfoil wake is described* by the ^ * 
boundary- layer-wake formulation. 

/ *" 

Wall Region: 
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Separation Bubble: Same as wall region except that 
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[/iny + /ivVl > 'Ll 

[V3y>/ \9x/ J = 6 


Airfoil Wake: Same as outer region 

Prandtl Number: 


Pr = 0.72 , Pr^ = 0.90 




The validity of. each of these models in the interaction and separated 
flow regions is highly suspect and requires verification by experiment. While 
it may be necessary to resort to more rigorous turbulent models in these 
regions, the simple models used here should permit some insight into the 
Influence of viscosity on such flow fields and are adequate for the early 
development stages of the computer codes. 
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Computational Time Step 


Six different computational time steps are used in the calculation; one 
for the Lx operator and one for the Ly operator in each of the three mesh 
regions. Each time step is determined by the CFL (Courant-Frledrichs-Lewy) 
and viscous stability requirements from the following relation: 

-c, . ... . ■ .. At - ~ • h ; 

I V| + a + (a/h) (v + e)/p “ 

where h is the appropriate mesh spacing, V is the appropriate velocity 
component, a is the local speed of sound, and a is a function of the mesh 
aspect ratio. In the fine mesh in the wake behind the airfoil, the eddy vis- 
cosity ,e is quite large, and the viscous stability criterion may govern the 
time step for , the Ly operator. To avoid. this undesirable restriction* and ' 
unneeded resolution of the wake, the entire fine mesh region' downstream of the " 
airfoil. is averaged and treated as part of the medium mesh. 

VALIDATION ■ ; .v' 

Because computations of this complexity have not previously been per- 
formed, the validity of the present code has been determined by comparison 
with established computations for certain specific regions of the flow field 
and. with experimental results obtained in the Ames High Reynolds Number 
Channel. 

' , Inviscid Flow Field 

The viscous terms were neglected and the inviscid flow was computed over 
a 6-percent circular-arc airfoil at a free-stream Mach number of 0.90. The 
surface pressure distribution is compared in figure 2 with a computation using 
the small-disturbance-theory program of E. M. Murman of Flow Research, Inc. 
(refs. 6 to 8 and unpublished information) and an Euler equation computation 
by R. VJ. MacCormack of Ames Research Center. Both of these computations 
employed special considerations at the shock which have not been incorporated > •:- 
in the present code. The agreement between the three computations (all of 
which solve the conservative form of the equations) is excellent, 

’ •• ■ r ^ SUL 


Attached Boundary Layer • 

The flow over an 18-percent-thick circular-arc airfoil at a free-stream 
Mach number of 0.775 and a chord Reynolds number of 2x10® was computed using 
eddy-viscosity model 1. The computed surface pressure distribution was input 
to the turbulent boundary-layer code of Marvin and Sheaffer of Ames Research 
Center (ref. 9 and additional information supplied by Marvin); which uses an 
implicit Crank-Nicholson method and a Cebeci-Smith eddy-viscosity model 
— (ref.; 10). Computed values of local skin friction from this boundary-layer, 
code are compared with the present calculation in figure 3. Ahead of the 
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separated-f low region the agreement . is excellent except for a small region 
near the leading edge where the Marvin-Sheaf f er code assumes a laminar to tur- 
bulent transition region. 

In addition to skin friction, boundary-layer velocity profiles computed 
by the present code are compared in figure 4 with the compressible form of the 
universal "law of the wall." Here the symbols represent the numerical solu- 
tion and are plotted at y"^ values corresponding to the mesh centers. All 
profiles are seen to have one point in 'the sublayer and adequately describe 
the log-law and wake-flow regions of the boundary layer. 

... Experimental '.Comparisons \ ^ , r ■ .v-,. ^ . 

Surface-pressure measurements, oil-flow studies, and shadowgraph studies 
were made, using an 18-percent-‘thick, circular-arc airf oil in the Ames High 
Re 3 rnolds Number Channel. Chord Rejmolds' numbers were varied between' 1 and 13' ■ 
million for free-stream Mach numbers between 0.711 and 0.788. In these ekper-'- 
iments, flow-field streamlines, determined from the present computer code, 
were used to design contoured tunnel walls, in an attempt to minimize tunnel 
Interference effects. 

Figure 5 shows a comparison between experimental and computed surface 
pressure distributions over the 18-percent-thick circular-arc airfoil at 
Moo “ 0*775 and Re^. = 2x10^. Three, computed distributions are shown in this 
figure as determined by the present code. One of these is an invisi'cid calcu- ' 
latlon, and the other two were determined using eddy-viscosity models 1 and' 2 '.^': 
All comparisons are in excellent agreement ahead of the interaction region.’ ' 
Both viscous computations show a marked improvement over the inviscld results 
in the interaction region. The solution with model, 2 shows the best agreement 
with experiment with respect to shock location and shock strength. Both vis- 
cous solutions indicate a shock- induced separation with the separation bubble 
extending into the wake. This is denoted in figure 5 by the pressure plateau 
downstream of the 80-percent-chord point. The experimental pressure distri- 
bution does not indicate this extensive separation region, and herein con- 
stitutes the greatest disagreement between computation and experiment. As 
will be shown later, at this particular choice of free-stream Mach number the 
experimental flow field can be highly unsteady, and .direct compa^rlspns in the. 
separated flow region may be invalid. • 

The first study performed using the code was to assess the influence of “ 
Reynolds number on the transonic flow field. Using the 18-percent circular 
arc at M<„ = 0.775, the chord Reynolds number was parametrically varied from 
1x10® to 6.67x10® and the flow. field computed using eddy-viscosity model 1. 

The results of this study are shown in figure 6 for the surface pressure dis- 
tribution. Included for reference is the inviscid solution. Three features 
are apparent in this figure. First the influence of Reynolds number on this 
flow field is small while the effect of viscosity is large. Second, as the 
Reynolds number is decreased, the shock strength decreases and the shock moves 
forward on the airfoil. And third, as the Reynolds number is decreased, the 
displacement effects in the separation region increase^ resulting in less pres- 
sure recovery near the trailing edge. The displacement effect ahead of the 
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Interaction region Is nearly the same for the entire Reynolds number range 
considered. 

A similar study was performed experimentally. Figure 7 shows the mea- 
sured surface pressure distributions over an 18-percent circular-arc airfoil 
at ^ = 0.750. Except for boundary- layer-transition effects at the lower 
Reynolds number, the experimental study bears out the same conclusions deter- , 
mined from the numerical results. In addition to Reynolds number effect, the 
Influence of free-stream Mach number was also experimentally studied. Fig- 
ure 8 shows measured surface pressure distribution over an 18-percent circular- 
arc airfoil at Re^. = lO^^lO^. Results for four free-stream Mach numbers are 
Included: Moo = 0.74, 0.76,- 0.77, and 0.79. At Moo = 0.74 -the flow is nearly 

subcrltlcal and the shock is very weak. For Moo ?= 0.76 the shock is much 
stronger but there is still a large pressure recovery at the trailing edge, .. 
Indicating small separation effects. For Moo = 0.77 there is a dramatic 
shift in pressure distribution over the aft portion of the airfoil. This - 
sudden jump is associated with the shock-induced separation merging with exis- 
tent trailing-edge separation, resulting in a large reverse flow region and . 
large boundary-layer displacement effects. . This phenomenon persists at the 
higher free-stream Mach numbers. This Mach number dependence of onset of mas- 
sive separation Is also a function of Reynolds number, which is discussed in • 
reference 11. 

To understand this discontinuous dependence on Mach member, shadowgraph 
movies were made of the flow over the aft portion of the airfoil as the Mach 
number was varied through the critical range. Photographs of selected frames 
of one such film are shown in figure 9 for an 18-percent circular-arc airfoil 
at Rej. = 7x10^. The Mach number variation is from 0.76 to 0.79 and results 
are shown for values of 0.76, 0.77, and 0.79. At Moo = 0.76 the flow is 
steady and there is some separation at the trailing edge. As the "Mach number 
is Increased, the flow becomes unsteady^ switching alternately from massive 
separation to fully attached or small trailing-edge separation. This 
unsteadiness is most probably an asymmetric phenomenon; hence it is not rea- 
sonable to expect to simulate this phenomenon with a symmetric, free boundary, 
code. When the Mach number reaches 0.79, the flow is again steady with mas- 
sive shock-induced separation. Mote that at the lower Mach number the shock 
is fairly weak and nearly normal to the airfoil surface. At the higher Mach 
number the shock is stronger and is definitely oblique — probably a lambda 
shock. 

To avoid comparisons in the unsteady flow regimes, subsequent computations 
were performed only for the lower and higher Mach numbers where the experi- 
mental flow is known to be steady. In addition, because of the superior per- 
formance indicated in figure 5 of eddy-viscosity model 2 versus eddy-viscosity 
model 1 in the interaction regime, model 1 has been dropped from further 
consideration. 

Oil-flow photographs of both the low and high Mach number experimental 
steady-flow regimes are shown in figure 10. At Moo = 0.76 the line at the 
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onset of trailing-edge separation is clearly indicated by the oil flow. The 
symmetry of the separation bubble is indicated in the right-hand photograph 
of the junction of the airfoil trailing edge and the tunnel side wall. At 
Moo = 0.79 the onset of shock-induced separation is indicated. The right-hand 
photograph of the junction of the airfoil trailing edge and the tunnel side 
wall clearly indicates the magnitude of the separation bubble. The photo- 
graphs in figure 10 also serve to indicate the two-dimensional and symmetric 
character of the flow. . 

For a chord Reynolds number of 4><10® the flow field over the 18-percent 
circular-arc airfoil was simulated using eddy-viscosity model 2 for = 0.742 
and Moo = 0.788, corresponding to the two extremes, of the experimental steady. . 
flows. ■ Figure 11 shows a comparison with experimental data of the surface 
pressure distributions for the two cases. In figure 11(a) the viscous solu- 
tion is an improvement over the Inviscid result but does not agree with the 
experiment in the shock— boundary-layer interaction region. The experiment 
indicates a stronger influence of viscosity than does the calculation, result- 
ing in a weaker, more smeared out shock. Near the trailing edge, however, the 
agreement is better, both distributions supporting the fact that there is only 
trailing-edge separation and both indicating the same level of pressure recov- 
ery. The reason for the substantial disparity in the interaction region is 
hot yet clear. It may be related to the fact that the shock is near the mid- 
chord and, for a chord Reynolds number of only 4^10^, the boundary layer may 
actually be transitional in that region. This could lead to the stronger 
vlscous-lnvlscld interaction effect indicated by the experiment. 

In figure 11(b) the viscous solution again is an improvement over the 
inviscid result. Here the agreement between experiment and viscous calcula- 
tion is good in the interaction region but is poor over the separation bubble. 
Unfortunately, there was no experimental data for Rej. = 4><10® at Moo = 0.788 
so data are shown for RSj. = 2x10® and 7x10®. It is expected that data for 
4x10® will fall within the envelope defined by these two limits. Note that 
both the experimental and calculated pressures tend to plateau over the sepa- 
ration bubble and indicate a similar extent over the aft portion of the airfoil. 
The computed pressures indicate larger recovery in this region and, as will be 
seen in the next figure, this is associated with the fact that the shock wave 
in the experiment is oblique while that simulated is nearly normal. The fun- 
damental reason' for this disparity is not yet clear but it is likely attribut- 
able to the inability of the simple eddy-viscosity model to support this com- 
plicated flow. 

Figure 12 shows the flow field detail over the aft portion of the airfoil 
for the two Mach numbers considered above. The top photographs are shadow- 
graphs of the experiment. For M^ = 0.742 the shock is weak and normal to 

the surface. Separation occurs only at the trailing edge and is small in 

extent. For M^g = 0.788 the shock is strong and oblique. Separation is 
shock Induced and extends into the wake. Immediately beneath the photographs 
are computed Mach number contours. For Moo = 0.742 the flow features are 

quite similar to the shadowgraph above. For Moo = 0.788 the shock is seen 
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to be nearly normal to the surface as opposed to being oblique in the shadow-, 
graph, but the magnitude of the separated flow region is similar to that shown 
in the shadowgraph. The bottom part of this figure contains computed velocity 
vector plots , showing the details of the separation bubbles and the effect.-.of 
the shock in retarding the flow. Subsequent experiments using a laser veloci- 
meter are planned to provide similar experimental data in this region. 


COMPUTER ‘requirements 


The present code was originally OTitten in FORTRAN and debugged on an 
IBM 360/67. using the interactive ..features , of the time sharing .system. Sub- 
sequently, it was run on the CDC 7600, requiring 2 to 10 hours per converged 
solution, depending on chord Reynolds number. Because of the long run times, 
the code was completely restructured and written in the vector-oriented ^es- 
developed CEU language. ^ This code, was translated to assembly language for 
the ILLIAC IV and to FORTRAN for the CDC 7600. The resulting FORTRAN code 
was further optimized, using COMPASS coded subroutines for aVl vector arith- 
metic operations. Resultant run times are now 0.8 to 4.2 hours per converged 
solution on the CDC 7600, and 0.6, to 3.0. hours on the ILLIAC IV. There is 
substantial room for speed Increase on the ILLIAC IV in that (1) the quoted 
run times were for the ILLIAC operating at 11.5 MHz Instead of the design 
speed of 15-16 MHz, (2) the ILLIAC was operated in non-overlap mode, and 
(3) each iteration was performed twice and the solutions compared before 
continuing with the computation. Each of these areas represents potential 
speed reductions of 0.72, 0.40, and 0.48, respectively, leading to an. overall 
potential speed reduction of 0.14. In this case, the present code would 
require from 0..08 to 0.42 hour, per , cony erged solution — an order of magnitude 
faster than the vector-coded CDC 7600. . , 


CONCLUDING, REMARKS 


In conclusion, a code has been written to simulate transonic turbulent 
flow fields over two-dimensional bodies of arbitrary configuration. At pres- 
ent, only algebraic eddy-viscosity models have been considered to achieve 
turbulence closure... With. .these -models the code yields valid solutions in the 
Invlscld floj^ field and in. the attached boundary .layer ahead of interaction 
regions. The validity of the numerical simulations in the shock— boundary- 
layer interaction region and in reverse flow regions is directly related to 
the . turbulence model. For the models considered thus far, comparisons with 
experiment have been less than good. This is to be expected since the models 
are developed from flat-plate incompressible boundary-layer data. Nevertheless 
the viscous solutions represent a considerable Improvement over inviscid solu- 
tions and do, predict the proper features of the flow. 


^Computational Fluid Dynamics (A FORTRAN-Based Language for the ILLIAC IV 
developed at Ames Research Center in 1973). 
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The code has been opt^ized for the CDC 7600 computer and requires from 

0.8 to 4.2 hours to simulate the transonic flow field over an 18-percent 
circular-arc airfoil. The code has also been written for the ILLIAC IV and : 
presently requires from 0.6 to 3.0 hours to simulate these transonic flows. 
These ILLIAC run .times have a potential reduction to 0.08 to 0.42 hour, if 
total’ advantage ' is taken of certain design features of the ILLIAC. 
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SYMMETRIC 

AIRFOIL 



Figure 1.— Computational control 
volume. 


Figure 2.— Invlscid pressure distribu- 
tion over a 6-percent circular-arc 
‘ airfoil; M<„ = 0.90. 



Figure 3.— Skin-friction distribution Figure 4.— Velocity profiles ahead of 
over an 18-percent circular-arc shock; = 0.775, Re = 2x10®. 

airfoil; ' = 0.775, Re^ = 2x10®. ^ 
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Figure 5.— Pressure; distribution . , 
over an 18-percent circular-arc 
airfoil; = 0.775, Re^ = 2x10^. 



0 


Figure 7.— Experimental pressure dis- 
tribution over an i8-percent 
circular-arc airfoil; Reynolds 
number effect, = 0.750. 



Figure 6.— Pressure distribution over 
an 18-percent circular-arc airfoil; 
Moo = 0.775. 



Figure 8.— Experimental pressure dis- 
tribution over an 18-percerit 
circular-arc airfoil; Mach number 

effect. Re = 10x10®. 
c 
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(bj Unsieady flo*^, oscillatory separation, A/« • 0. 77 



(c) Steady florv , shock-induced separation , 79 


Figure 9.— Photographs of boundary-layer separation from a shadowgraph movie; 

t/c = 0.18, Re^ = 7x10^. 
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